#include <stdio.h>
#include <assert.h>
#include <sys/time.h>
#include "xml_config.h"
#include <gsl/gsl_math.h>
#include <gsl/gsl_const_mksa.h>

// Open MultiProcessing
#ifndef SINGLE_THREAD //Declare this to disable multithreaded behaviour
#include <omp.h>
#endif

#include "cdecl.h"
#include "func_EHw.h"


/**
 * This file contains the glue-up program of the whole.
 */


/** 
 * Reports how much time have passed since programm start and how 
 * much time it is at themoment of calling
 * @param start - moment, since what time difference is calculated
 * @param message - text to identify moment (will be printed to stderr)
 */
static void report_time(const struct timeval * start, const char * message){
	struct timeval now;
	gettimeofday(&now, NULL);
	
	
	time_t dsec = now.tv_sec - start -> tv_sec;
	suseconds_t dusec; //It's unsigned, so we can't let it go negative
	if(now.tv_usec < start -> tv_usec){
		dsec -= 1;
		dusec = 1000000 - (start -> tv_usec - now.tv_usec);
	}else{
		dusec = now.tv_usec - start -> tv_usec;
	}
	fprintf(stderr, ">>%.6fs - %s\n",(double)dsec + 1e-6*(double)dusec, message);	//It's OK to cast to int, since we won't overflow it in practical situations ^^
}

/**
 * Simple random generator with internal state in 'state' variable.
 * Generates a value from [0,m] diapason.
 * @param state - generator internal state, shall be set to something 
 * before generation.
 * @param m - generate numbers up to m inclusive
 * @return low quality pseudorandom number
 */
static uint32_t myrandom(uint64_t * state, uint32_t m){
	// Deciding about mask
	uint32_t imask = ~0;
	do{
		imask <<= 1;// Shift mask 1 bit left
	}while( ~imask < m );
	uint32_t mask = ~imask;
	// Now we have a minimal mask that covers all non-zero bits of m
	uint32_t ans;
	do{
		* state = * state * 8216778763LL - 2148773723LL;
		 ans = (* state >> 20)  % UINT32_MAX;
		 // Applying mask
		 ans &= mask;
	}while(ans > m);// After some tries we'll get it
	return ans;
}


/**
 * Generates a random shuffle using Durstenfeld aka Fisher-Yates 
 * algorythm. Numbers from {0,...,N-1} are randomly distributed in
 * N cells of shuffle array.
 * The shuffle is supposed to be used for OpenMP spectrum
 * calculation planing between concurrent threads. Without such planning
 * some thhreads tend do do the work faster than the other.
 * @param N - length of the shuffle
 * @return pointer to malloc-ed buffer with shuffle data or NULL, if
 * memory could not be allocated.
 */
static uint32_t * shuffle(uint32_t N){
	// Random seed based on UNIX time
	uint64_t random_state = time(NULL);// Unix time to seed
	uint32_t * values = malloc(sizeof(uint32_t)*N);
	if(NULL == values){
		return NULL;
	}
	// Alloc OK
	uint32_t ksi;
	for(ksi = 0; ksi < N; ksi += 1){
		values[ksi] = ksi;
	}
	for( ksi = N-1; ksi > 0; ksi -= 1 ){
		/* for all positions from max down to second from beginning
		 * inclusive
		 **/
		/* ksi is an index, above which elements of permutation have
		 * allready been chosen and will be not modified
		 **/
		uint32_t choice = myrandom(&random_state, ksi);
		uint32_t e1, e2; /* Things to exchange */
		e1 = choice;
		e2 = ksi;
		{/* exchange e1 and e2 */
			uint32_t tmp;
			tmp = values[e1];
			values[e1] = values[e2];
			values[e2] = tmp;
		}
	}
	return values;
}


struct clean_up_struct {// Objects that need to be fereed properly, NULL disables freeing
	struct configuration * cfg;
	complex * signal_w;
	struct discharge * reflections;
};

static void clean_up (struct clean_up_struct * cl){
	if(NULL != cl -> cfg){
		configuration_free(cl -> cfg);
	}
	free(cl -> signal_w);// NULL OK
	trash_discharge_list( cl -> reflections );
}

#define PRECISION 0.001
#define T_CUTOFF_PRECISION 0.01
#define W_CUTOFF_PRECISION 0.001 //0.002 used to be ok
#define OMEGA_MAX_RESERVE 2.0 //Make signal max spectrum Omega_max this times more than obtained from estimations
#define T_RESERVE 1.7 //Same with max signal duration T
 // 1.5 1.5 used to be good here
 

int main(int argc, char ** argv){
	//Remembering start time
	struct timeval start;
	gettimeofday(&start, NULL);
	report_time(&start, "Start");
	
	//Preapring clean-up
	struct clean_up_struct cl;
	//Init all to NULL until set to something
	cl . cfg = NULL;
	cl.signal_w = NULL;
	cl.reflections = NULL;
	// Init all to NULL here
	
	
	report_time(&start, "reading config");
	// Reading configuration from file
	struct configuration cfg;
	int xml_read_status = xml_config_parser(&cfg, 0); // Reading from stdin
	if(xml_read_status){
		fprintf(stderr, "Failed to read XML config. Parser returned %i\n", xml_read_status);
		clean_up( & cl );
		return 1;
	}
	double eps = PRECISION;
	double cweps = W_CUTOFF_PRECISION;
	double cteps = T_CUTOFF_PRECISION;
	cl.cfg = & cfg; // For clean-up
	
	//Creatng mirrored discharges if that's necessary
	struct discharge * reflections = NULL; // yet empty
	if(cfg.has_mirror){
		report_time(&start,"creating discharge reflections");
		struct discharge * tmp_root = NULL; // temporary reflections list root
		struct discharge * dis;
		for(dis = cfg.discharges; dis != NULL; dis = dis -> next){
			struct discharge * refl = malloc(sizeof(struct discharge));
			if(NULL == refl){
				//Wooops!
				fprintf(stderr, "Out of memory when making discharge reflections (1)\n");
				trash_discharge_list(tmp_root);
				clean_up( & cl );
				return 5;
			}			
			int mres = discharge_mirror(dis, refl);
			if(mres){
				fprintf(stderr, "Out of memory when making discharge reflections (2)\n");
				trash_discharge_list(tmp_root);
				clean_up( & cl );
				free(refl);
				return 25;
			}
			//OK mirroring, linking into the list
			refl -> next = tmp_root;
			tmp_root = refl;
		}
		// Publishing tmp_root
		reflections = tmp_root;
		cl.reflections = reflections; // For clean-up
	}
	
	// Let's make array of pointers to all available discharges for convinience
	int nd = 0;
	{
		// How many do we have?
		struct discharge * dis;
		for(dis = cfg.discharges; dis != NULL; dis = dis -> next){
			nd += 1;
			assert(nd > 0);
		}
		for(dis = reflections; dis != NULL; dis = dis -> next){
			nd += 1;
			assert(nd > 0);
		}
	}
	struct discharge * discharge_pointers[nd];
	{
		//Initializing the pointers
		int ksi = 0;
		struct discharge * dis;
		for(dis = cfg.discharges; dis != NULL; dis = dis -> next) {
			discharge_pointers[ksi] = dis;
			ksi += 1;
		}
		for(dis = reflections; dis != NULL; dis = dis -> next) {
			discharge_pointers[ksi] = dis;
			ksi += 1;
		}
	}
	
	
	report_time(&start,"deciding about fft parameters");
	// Deciding about fft parameters
	double Omega_max, T;
	{
		//Maximum Omega, geneally, deends on trajectory also. 
		//Let's estimate it, however, by current spectrum taking 
		//trajectory in consideration not.
		
		// With max duration it shall be all OK
		struct discharge * dis;
		double maxw = 0, maxt=0;
		
		int ksi;
		for(ksi = 0; ksi < nd; ksi += 1){
			struct discharge * dis = discharge_pointers[ksi]; // All here!
			double thismaxw = current_get_max_w( & dis -> crnt, cweps);
			if(thismaxw > maxw){
				maxw = thismaxw;
			}
			double L = trajectory_len_calc(& dis -> trlen, 1.0); // Whole trajectory length
			double thismaxt = dis -> t_start + 
					current_get_duration(& dis -> crnt, L, cteps) +  // Current delay
					trajectory_max_dist(& dis -> trj, cfg.r) / GSL_CONST_MKSA_SPEED_OF_LIGHT;//propagation delay (we are in ~ air)
			if(thismaxt > maxt){
				maxt = thismaxt;
			}
		}
		Omega_max = maxw * OMEGA_MAX_RESERVE;
		T = maxt * T_RESERVE;
	}
	double Omega = 2*M_PI / T; // This is a different value from Omega_max, but looks similar.
	
	// fft size
	uint32_t N;
	{
		double NN = 2*Omega_max/Omega;
		if(NN > UINT32_MAX){
			fprintf(stderr, "fft dimension too big (futuristic)\n");
			clean_up( & cl );
			return 2;
		}
		N = (uint32_t)NN;
		assert(N > 60);//this is needed for linear addition removal, at least
		if (! (N % 2)){
			N -= 1;// We want to have odd number of knots
		}
	}
	double h_t = T / N; //Time step of signal after fft 
	
	report_time(&start,"calculating spectrum");
	// Allocating place to store spectrum in 
	complex * signal_w = malloc( N * sizeof(complex)); // signal spectrum
	if(NULL == signal_w){
		fprintf(stderr, "Out of memory when allocating space for fft coefficients :(\n");
		clean_up( & cl );
		return 3;
	}
	cl.signal_w = signal_w; // Remembering for clean-up
	
	
	// The sweet for all frequencies cycle 
	uint32_t task_id, mu = N / 2;// fft vector index k and positive -> negative spectrum cut index mu
	{
		// Preparing a random shuffle
		uint32_t * task_shuffle = shuffle(mu + 1);
		if(NULL == task_shuffle){
			fprintf(stderr, "Out of memory when allocating task shuffle :(\n");
			clean_up( & cl );
			return 31;
		}
		
		uint32_t N_failed = 0;//Failed to calculate integral for so many frequencies 
		#ifndef SINGLE_THREAD
		#pragma omp parallel for reduction(+:N_failed) 
		#endif
		for(task_id = 0; task_id <= mu; task_id += 1){
			uint32_t k = task_shuffle[task_id];
			if(cfg.field_class == FIELD_E && k == 0){ 
				// For E field we have singularity in w = 0. The teory says,
				// that if we blank the k=0 coefficient like here and 
				// calculate signal without it, we will get (signal(t) + C*t + D)
				// as a result. The unwanted C*t + D will be removed later
				signal_w[0] = 0;
				continue;
			}
			double w = Omega * k;
			
			signal_w[k] = 0;
			int ksi;
			for(ksi = 0; ksi < nd; ksi += 1){// For all discharges
				struct discharge * dis = discharge_pointers[ksi];
				complex eh_res;
				int eh_status = func_EHw_calc(cfg.field_class, eps, cfg.r, cfg.n, w, & dis -> trj, & dis -> trlen, & dis -> crnt, & eh_res);
				
				if(eh_status){
					fprintf(stderr, "Failed calculating spectrum @ w = %f\n", w);
					N_failed += 1;
					signal_w[k] = GSL_NAN + GSL_NAN*I;
					break;
				}
				
				signal_w[k] += eh_res * cexp(I*w*(dis -> t_start)); // Considering delayed discharge start
				// Report progress
				if(!(k % 20)){//not to flood tty
					const char * format = "%li: f=%.2f Hz\t(thread %i)";
					int tid;
#ifdef SINGLE_THREAD
					tid = 0;// Allways zero
#else
					tid = omp_get_thread_num();
#endif
					double f = w / (2 * M_PI);
					int ls = snprintf(NULL, 0, format, k, f, tid);
					char msg[ls + 1];
					snprintf(msg, ls + 1, format, k, f, tid);
					report_time(&start,msg);
				}
			}// end for all discharges
			
			// Attenuating spectrum for use with fft
			signal_w[k] *= Omega;
		}// Done positive specrtm part
		
		free(task_shuffle); // Deallocating before next return
		
		//Checking If we have had errors
		if(N_failed){
			fprintf(stderr, "Spectrum calculation failed %ui times\n", N_failed);
			clean_up( & cl );
			return 4;
		}
	}
	
	{
		uint32_t k;
		// Using real signal spectrum feature: S(-w) = conj(S(w))
		for(k=mu+1; k < N; k += 1){
			uint32_t k_m = N - k;
			signal_w[k] = conj(signal_w[k_m]);
		}
	}
	
	//DEBUG!
#ifdef DEBUG
	{
		FILE * f;
		f=fopen("sd.dat", "w");
		uint32_t o;
		for(o=0;o<N;o+=1){
			fprintf(f, "%e, %e\n", creal(signal_w[o]), cimag(signal_w[o]));
		}
		fclose(f);
	}
#endif
	
	
	report_time(&start,"doing fft");
	//Applying FFT
	int fft_res = numba_cruncha_fft(signal_w, N);// now signal_w holds time-domain values
	if(fft_res){
		fprintf(stderr, "Out of memory when calculating fft :(\n");
		clean_up( & cl );
		return 7;
	}
	
	complex * signal_t = signal_w; // Just copying pointer for convinience
	// Correcting k*t + b if that was E field
	if(cfg.field_class == FIELD_E){
		report_time(&start,"Removing a random k*t+b functional addition from the signal.");
		long takelast=10;
		assert(2*N > takelast);/* to keep of out of bounds event */
		double ff[takelast];
		double tt[takelast];
		long j;
		for(j=N-takelast;j<N;j++){
			ff[j-N+takelast]=creal(signal_t[j]);
			tt[j-N+takelast]=j*h_t;
		}
		double k,b;
		linavg(tt,ff,takelast,&k,&b);/* calculating k*t+b approximation */
		fprintf(stderr,"k=%g,b=%g\n",k,b);
		double E_0_=creal(signal_t[0]); /* this will be substracted 
		* from E(t), so that E starts from zero (just a beauty) */
		for(j=0;j<N;j+=1){/* adding -k*t-E_0_ function to field */
			double t=h_t*j;/* current time */
			signal_t[j]=signal_t[j]-k*t-E_0_;
		}
	}
	
	
	report_time(&start,"exporting signal");
	// Exporting signal
	printf("%% time,\tsignal\n");
	{
		uint32_t k;
		for(k = 0; k < N; k+=1){
			double S_t = creal(signal_t[k]);
			double t = k * h_t;
			printf("%e, %e\n", t, S_t);
		}
	}
	
	
	// Cleaning up
	clean_up( & cl );
	
	
	// Quk-ik
	report_time(&start,"Done.");
	fprintf(stderr,"Quk-ik, Ja Kuigorozsh.\n");
	fprintf(stderr,"     ^  ^            \n");
	fprintf(stderr,"     \\00//=====-      \n");
	fprintf(stderr,"   //\\\\/ \\ ==-        \n");
	fprintf(stderr,"   // \\   \\           \n");
	fprintf(stderr," //     \\ \\   //\\/\\/\\/\n");
	fprintf(stderr,"         \\ \\/ /       \n");
	fprintf(stderr,"          ----        \n");
	fprintf(stderr,"Davajte mne rabotu, a neto\n");
	fprintf(stderr,"vse porushu-polomaju...\n");
	
	return 0;
}
